Screening in Ionic Systems: Simulations for the Lebowitz Length 
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Simulations of the Lebowitz length, £l(T, p), are reported for the restricted primitive model hard- 
core (diameter a) 1:1 electrolyte for densities p < 4p c and T c < T <■ 40T C . Finite-size effects 
are elucidated for the charge fluctuations in various subdomains that serve to evaluate £l. On 
extrapolation to the bulk limit for T > 10T C the low-density expansions (Bekiranov and Fisher, 
1998) are seen to fail badly when p > yjjpc (with p c o? ~ 0.08). At higher densities £l rises above 
the Debye length, £d oc v/T/ p, by 10-30% (upto p ~ 1.3p c ); the variation is portrayed fairly well by 
generalized Debye-Hiickel theory (Lee and Fisher, 1996). On approaching criticality at fixed p or 
fixed T, £,l(T, p) remains finite with ~ 0.30a ~ 1.3£fj but displays a weak entropy-like singularity. 
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Understanding the thermodynamic and correlation 
properties of ionic fluids has challenged both theory and 
experiment Q. Typical electrolytes exhibit phase sep- 
aration that is analogous to the gas-liquid transition in 
simple fluids, albeit at rather low temperatures when ap- 
propriately normalized. However, the long range of the 
Coulomb interactions has hampered understanding espe- 
cially near criticality Q. One crucial aspect is Debye- 
Hiickel screening. For a d-dimensional classical fluid 
system with short-range ion-ion potentials beyond the 
Coulomb coupling z a z T q 2 /r d ~ 2 (where z a is the valence 
of ions of species a and mole fraction x a while q is an ele- 
mentary charge) , the charge-charge correlation function, 
Gzz(r;T,p), decays as exp[-\r\/£ z ,oo(T, p)] (see, e.g., 0, 
y]): the asymptotic screening length, S,z.oo, approaches 
the Debye length £d = (k^T / Airz 2 q 2 p) 1 / 2 when the over- 
all ion density p approaches zero (with z\ — ^2 z 2 x a 

Hi). 

By contrast, at a critical point of fluid phase sepa- 
ration, the density-density (or composition) correlation 
length, £,n,oo{T, p), diverges, as do all the moments of 
Gnn{i",T, p). What then happens to charge screening 
near criticality? This question was first posed over a 
decade ago [4| and has been addressed recently via the 
exact solution of (d > 2)-dimensional ionic spherical mod- 
els As anticipated [4(b)], the issue of ± ion symme- 
try proves central. However, spherical models for fluids 
display several artificial features (e.g., infinite compress- 
ibilities on the phase boundary below T c ; parabolic co- 
existence curves, (3 = i; etc.). Accordingly, understand- 
ing screening near criticality for more realistic models 
remains a significant task. 

To that end we report here on a Monte Carlo study 
of the restricted primitive model (RPM), namely, hard 
spheres of diameter a carrying charges q± = ±q (so 
that z + = —z- = 1, x+ = X- = i). Grand canon- 
ical simulations have been used and, to accelerate the 
computations, a finely discretized (£ = 5 level) lattice ver- 
sion of the RPM has been adopted 0- For this system 
the critical behavior is well established as of Ising-type 



with T* = k B T c a/q 2 ~ 0.05069 and p* ee p c o? ~ 0.079 0. 
Furthermore, it has been demonstrated that for £ > 3 
the fine-lattice discretization does not qualitatively affect 
thermodynamic or finite-size properties Q. 

Ideally one would like to calculate £n,oo(T,p) and 
£,z,oa(T, p) near criticality; but even in nonionic model 
fluids, obtaining £at j0 o via simulations is hardly fea- 
sible. Nevertheless, the low-order moments Mj^fc = 
J \r\ k GNN (r)d d r for k = 0, 1, 2, • • • , are accessible and, 
by scaling, all the = {M N ^ k /M N ^) x / k for k > 

diverge like £jv,oo- However, for charges the Stillinger- 
Lovett sum rules 0,0 dictate Mz.o = (so that Gzz(f) 
is not of uniform sign) while the second moment satisfies 
Mz,2 = 2z|<7 2 p£q which is fully analytic through (T c , p c ). 
On the other hand, the first moment of Gzz{r) is known 
to be intimately related to charge screening via the 
so-called "area law" of charge fluctuations. 

To explain this, consider a regular subdomain A with 
surface area and volume | A| , embedded in a larger do- 
main, specifically say, the cubical L d simulation box. If 
Qa is the total fluctuating charge in A, electroneutrality 
implies (Qa) = 0; but the mean square fluctuation, (Q'a), 



will grow when |A| increases 

12 



In the absence of screen- 



A ) is asymptotically 
. This was first ob- 



ing one expects (Qa) ~ |A|; however, in a fully screened, 
bulk (L — > oo) conducting fluid 
proportional to the surface area 
served by van Beijeren and Felderhof and later proven 
rigorously by Martin and Yalcin [8| . Following Lebowitz 
[8| one may then define a screening distance proportional 
to Mz,\{T, p), which we call the Lebowitz length, p) 
2\ via' 



(Qi)/AA^C dP Z 2 2 q 2 ^(T,p) 



|A| 



(1) 



Note 



where Cd is a numerical constant with cj, = \- 
that, since Gzz{ r ) is n °t necessarily of uniform sign, 
£,h(T,p) oc Mz,i(T,p) might diverge at T c even though 
the second moment Mz,i oc remains finite! 

Clearly, by simulating (Q\) in various subdomains one 
may, as we show here, hope to calculate the Lebowitz 
length. To our knowledge no numerical results have been 
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reported previously for d — 3 although Levesque et al. 
y| presented a study (above criticality) for d = 2. An 
exact low density expansion proves that £l/£d — * 1 
when p — > and corrections of order P 1 ^ 2 , phip and 
p have been evaluated. This analysis [2| also served to 
validate the generalized Debye-Hiickel (GDH) theory for 
the correlations 01 f° r small p. 

The GDH theory, however, did not generate a p\np 
term: nevertheless, as we find here, the exact expan- 
sion fails at very low densities — around p c /lQ even for 
T ~ 10T C — while GDH theory provides a reasonable 
estimate of £,h(T,p) at higher densities: see Fig.|21below. 
Furthermore, our calculations show that £l remains finite 
at criticality, exceeding by only 33%. Nonetheless, 
the Lebowitz length does exhibit weak singular behavior 
that, in accord with general theory, matches that of the 
entropy. 

The first serious computational task is to understand 
the finite-size effects resulting from the L x L x L sim- 
ulation box with periodic boundary conditions. Each 
simulation at a given (T* , p* ) yields a histogram of 
the total fluctuating charge Q\ for 24 different subdo- 
mains A. We have used: six small cubes of edges XL 
with A = 0.3, 0.4, •••,0.8; seven 'rods' of dimensions 
XLxXLx L with A = 0.2, . . . , 0.8, four 'slabs' of di- 
mensions XL x L x L with A = 0.2, • • • , 0.5; and seven 
spheres of radius R = XL with A = 0.15-0.45 in in- 
crements AA = 0.05. To minimize correlations between 
these various subdomains, they have been located as far 
apart as feasible. 

While the area law for the charge fluctuation, (Q\), is 
rigorously true for L — * oo followed by A — > oo, it is by no 
means clear how it will be distorted for a finite subdomain 
A embedded in a finite system. To understand this Fig.Q] 
presents {Q\), normalized by q 2 , for the six cubic subdo- 
mains as a function of the reduced area Aa/L 2 at selected 
temperatures and densities for box sizes L* = L/a = 6 
and 12. Surprisingly, at high temperature and moderate 
density (T* = 0.5 ~ 10T*,p* = 0.08 ~ p* c ), the area law 
is well satisfied for A < 0.7 even for small systems. For 
L* = 6 the data point for A = 0.8 deviates strongly from 
the linear fit (dashed line) owing to finite-size effects: in- 
deed, electroneutrality dictates that (Q A ) should vanish 
when A — > 1, corresponding to Aa/L 2 = 6. At low densi- 
ties around |p c , the Debye length £d oc \jTj p becomes 
large but nevertheless we see that the area law is still 
well satisfied. Furthermore, the area law is found to hold 
even near criticality: see the lowest plot. Note, however, 
that the linear fits to the data do not pass through the 
origin. This reflects finite-size effects which are discussed 
further below. 

Combining JQl with the observations illustrated in 
Fig. Q we conclude that charge fluctuations in the cu- 
bic subdomains are well described by 

(Q 2 A (T, p; L)) = Aq(T, p; L) + ± W 2 £ L (T, p; L)A A , (2) 
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FIG. 1: Reduced charge fluctuations, (Q\) /q 2 for given (T, p) 
in cubes A of edges XL vs. reduced area, Aa/L 2 = 6A 2 . 

where the intercept Ao(T, p; L) need not vanish. The (fit- 
ted) linear slope serves to define the finite-size Lebowitz 
length, ^l{T, p; L), which should approach the bulk value, 
£l (T, p) . But by what route? 




FIG. 2: Quadratic fits to finite-size Lebowitz length data for 
sizes up to L* — 24 at T* = 0.5 and various densities. 

To answer this question consider Fig. [3 which displays 
£l(T,p;L) vs. 1/L* for T* — 0.5 at various densities. It 
is rather clear that £l(T, p; L) approaches its bulk limit 
as 1/L. This can be understood by recalling the Lebowitz 
picture in which the uncompensated charge fluctua- 
tions in a subdomain arise only from shells of area A\ 
and thickness of order £ L • By invoking the screening of 
Gzz{r) one can see that A£l =£l(£) — £l(oo) for smooth 
subdomains decays as \j_L 2 ■ Indeed, by this route van 
Beijeren and Felderhof [8j showed explicitly that fluctu- 
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ations in a sphere of radius R (in an infinite system) ap- 
proach their limiting behavior as 1/R 2 . For spheres in fi- 
nite systems, we observe similarly that £l(£) approaches 
the bulk value as 1/L 2 . However, for cubes — which have 
edges and corners — and rods with edges, £l(L) gains a 
lower order, 1/L term as seen in Fig. EI (The intercept 
Aq(L) in (0) is, correspondingly, found to vary as L.) 
On the other hand, for slabs, lacking edges and corners, 
we find that obtained via Q approaches the limit 

exponentially fast. 



1.25 




FIG. 3: Density variation of the bulk Lebowitz length extrap- 
olated from various subdomains at T* = 0.5. The dashed, 
solid and dotted plots represent GDH theory and ap- 

proximants exact at low density: see text. 

Having established the finite-size behavior, let us ex- 
amine £l(T,p) on the T*=0.5 isotherm, well above 
T c . Figure |2l shows estimates extrapolated from cubes, 
spheres and slabs. At moderate densities systems up to 
L* = 16 suffice but for p* < 0.025 we went up to L* = 24. 
The results may be compared with GDH theory 
(dashed curve) and approximants which reproduce the 
exact low-density expansion known to order p ,2]. For 
the latter we adopt 



1.0] 



&(T,p)[l + oi(T)p* 

€D(r,p)/[i-oi(r)p' 



oa(T)p*lnp*], (3) 
-<z 2 (T)p* In p*j, (4) 



shown in Fig. [31 as solid and dotted curves, respectively, 
where a\(T) and 02 (T 1 ) follow from Q. The simulations 
agree well with the low-density expansion but only up 
to p* ~ 0.005; thereafter £l rises above the Debye length 
much more slowly. By contrast, GDH theory captures the 
overall behavior of £l{T, p) over a broad density range, 
representing the numerical estimates to within a few per- 
cent at moderate densities, 0.01 < p* < 0.10, where no ex- 
act results are available. 

In the critical region the first question is the finiteness 
of £l(Ic, p c ). To answer we study £l on the critical iso- 



chore p = p c as T^T C . Figure l4l \Ti\ reveals that Cl/o 
falls increasingly rapidly when T* drops from ~ 0.5 but 
clearly attains a finite nonzero value at T c that exceeds 
f£/a~ 0.2260 @. Owing to the relatively strong finite- 
size dependence of £l and the excessively large computa- 
tional requirements near (T c ,p c ), reliable extrapolation 
to L = oo is difficult. Nevertheless we may test for the 
nonanalytic behavior expected in any finite quantity |l2| . 
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FIG. 4: Lebowitz length for L* = 8 and 12 on the critical 
isochore compared with the Debye length. 

On general grounds |l2l | weak, entropy-like behavior is 
predicted. Thus temperature derivatives at p = p c should 
diverge like the specific heat, namely as 



pC v /k B ~ A+/t a + A° 



(5) 



when i=(T - rj/T c 



,3 



I) 



0.50 ±0.07 13] with, 



where a ~ 0.109 and 
via a rough fit, A a 3 ~ 



—0.37. A direct comparison for finite L of 9(£l/£d)/<9T 
with the specific heat is shown in Fig. El ^3 - Bearing in 
mind the lack of £l data near T c and its imprecision, the 
resemblance of the two plots is striking: we accept it as 
confirmation of the anticipated singularity. 

Complementary nonanalytic behavior should arise on 
the critical isotherm as the reduced chemical potential 
p* = [p-p, (T)]/kBT varies. This is borne out by the 
plots in Fig. El of d(£ h /&)/dp* and (d(p*U*)/dp*)/ p* k 
with k = i , where U* (T, p) is the configurational cn- 
e po 

m. 

by scaling, diverge as 1/ 



2' 

ergy per particle; the power p* K represents a conve- 
nient "fc-locus factor" |l5|. In the bulk limit both func- 
tions should, by scaling, diverge as 1/1 u — a c \^ with 



^ = (l-/3)/09 + 7)^0.43 
Returning to the isochore p 



= p c , theory indicates 



9* 1 



e 2 t 2 + 



^{T)=^[l + e a t 1 - a + e 1 t + , 

where 8 ~ 0.52 is the leading correction exponent 0| . By 
making allowance for the L-dependence and fitting over 
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FIG. 5: (a) Temperature derivative of reduced Lebowitz 
lengths and (b) specific heats on the critical isochore. The 
dashed curve approximates the bulk specific heat |13| . 
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FIG. 6: Derivatives on the critical isotherm of (a) the reduced 
Lebowitz lengths and (b) the energy densities with respect to 
the chemical potential fj* where 1. 36218: see text. 



various ranges above T c we conclude ££ ~ 0.30a and, with 
less confidence, e a ~ 2.6 ± 0.2 and e\ ~ — 2.2 ± 0.3. 

In summary, the Lebowitz screening length, £,l(T, p), 
has been studied for the restricted primitive model elec- 
trolyte via grand canonical Monte Carlo simulations of 
the charge fluctuations in subdomains. The correspond- 
ing area law that is asymptotically valid for large sub- 
domains Q holds surprisingly well even in small simula- 
tion boxes, L < I2a. Finite-size effects can be understood 
so that the bulk, L — > oo limit may be extracted by ex- 
trapolation vs. 1/L for cubic subdomains and l/L 2 for 
spheres while the effective, finite-size Lebowitz lengths 
for slabs converge exponentially fast. Evaluation of £l 
for T > 10T C over densities from 0.03p c to Ap c reveals that 



the exact low-density expansions J2J are effective only for 
plSjoPc whereas GDH theory [nj reproduces well the 
general trends. Finally, £l remains finite at criticality but 
exhibits weak, entropy-like singularities on approaching 
(T c , p c ). This is the first time that charge-charge corre- 
lations and a strongly state-dependent screening length 
have been studied by simulations close to criticality. 
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